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This paper focuses on the problem of growing multiplex networks. Currently, the results on the 
joint degree distribution of growing multiplex networks present in the literature pertain to the case of 
two layers, and are confined to the special case of homogeneous growth, and are limited to the state 
state (that is, the limit of infinite size). In the present paper, we obtain closed-form solutions for 
the joint degree distribution of heterogeneously growing multiplex networks with arbitrary number 
of layers in the steady state. Heterogeneous growth means that each incoming node establishes 
different numbers of links in different layers. We consider both uniform and preferential growth. We 
then extend the analysis of the uniform growth mechanism to arbitrary times. We obtain a closed- 
form solution for the time-dependent joint degree distribution of a growing multiplex network with 
arbitrary initial conditions. Throughout, theoretical findings are corroborated with Monte Carlo 
simulations. The results shed light on the effects of the initial network on the transient dynamics of 
growing multiplex networks, and takes a step towards characterizing the temporal variations of the 
connectivity of growing multiplex networks, as well as predicting their future structural properties. 


I. INTRODUCTION 

The network framework is widely used for studying 
complex systems and their properties, through mapping 
units onto nodes and interactions onto links. The con¬ 
ventional conceptualization consists of a single graph, 
while many real-world systems exhibit different types 
of interactions between constituents. Such networked 
systems in which units have heterogeneous types of in¬ 
teractions can be mathematically modeled under the 
framework of multiplex networks. In these settings, 
units (nodes) are members of distinct networks simul¬ 
taneously. For example, it would be an oversimplified 
depiction to talk about the social network of a given 
group of people, because the same pair of individuals can 
have distinct types of relations at the same time: they 
can be kins, friends, coworkers, they can be connected on 
the social media, etc. In other words, there can be mul¬ 
tiple networks between the same set of people. This is 
the basic rationale behind the multiplex representation 
of networked systems. 

The multiplex framework envisages different layers 
housing different types of links between the same set 
of nodes (so there is one set of nodes, multiples sets 
of links). For example, we can take a sample of indi¬ 
viduals and constitute a social media layer, in which 
links represent interactions on social media, a kinship 
layer, a geographical proximity layer, and so on. Exam¬ 
ples of studies that have conceptualized real networked 
systems under the multiplex framework include citation 
networks mm, online social media mu, interbank net¬ 
works [3, airline networks [5], scientific collaboration 
networks mm, and web of connections and interactions 
in online games [8j. 

Theoretical tools for characterizing and quantifying 
the properties of multiplex networks are generalizations 
of the single-layer scenario to multiple layers pin]. For 
example, the adjacency matrix is generalized to the ad¬ 


jacency tensor, whose ijk element is the weight of the 
link from node i to node j in layer k. Similarly, all 
nodal attributes which were scalars in the single-layer 
picture (such as various types of centralities, degree, 
and clustering) are generalized to vectors in the mul¬ 
tiplex scenario. These new theoretical measures enable 
studying various phenomena on top of multiplex net¬ 
works analytically. Examples include epidemics [I1II31, 
pathogen-awareness interplay [Til [T5] . percolation pro¬ 
cesses nadiiiiTi, random walks [HilS], evolution of 
cooperation [2(11123] . diffusion processes [24] and social 
contagion Eor thorough reviews, see Emanz]. 

Since many real networks are growing in size, growing 
multiplex networks have also attracted attention in the 
literature. The mean-field approach is a potent method 
for investigating the temporal evolution of the degrees of 
individual nodes, extracting its asymptotic behavior in 
order to find the asymptotic (tail behavior) degree dis¬ 
tribution of each of the individual layers. This approach 
is undertaken in [2HHSI]- An alternative approach of 
tackling network growth problems is the rate equation 
approach, undertaken in [2H1 IMl E2 • The rate equation 
enables solving for the joint degree distribution of the 
system, so that we can obtain the fraction of nodes that 
have a given degree vector across layers. 

Previous results on the joint degree distribution of 
growing multiplex networks—attainable through the 
rate equation approach—are confined to the case of ho¬ 
mogeneous two-layer growth, where the number of links 
established by each new incoming node is the same across 
layers [28l|29]. Note that the possibility of heterogeneous 
growth is envisaged in [2H| (in the Supplemental Mate¬ 
rial therein), and their implications for the mean-field 
scenario are correctly alluded to. In the present paper, 
we consider different rates of link growth across layers ex¬ 
plicitly, and obtain the joint degree distributions. More¬ 
over, we extend the problem to general M layers. To our 
knowledge, no solution for the joint degree distribution 
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of growing multiplex networks with arbitrary number of 
layers exists in the literature. Furthermore, the existing 
results on the joint degree distribution are confined to 
the steady-state, that is, the limit as t —)■ oo where the 
network has infinite size. The evolution of the joint de¬ 
gree distribution at arbitrary times is hitherto unknown. 
In this paper, we take a step towards alleviating these 
two gaps. 

The assumption of link growth heterogeneity is moti¬ 
vated by the empirical studies on the structure of multi¬ 
layered interacting systems, which report that different 
layers generally exhibit different connectivity patterns 
(consequently, different average degrees, and other struc¬ 
tural properties). For example, in [5], the connectiv¬ 
ity structure of the players of a massive online game is 
mapped onto a multiplex network of six layers, and the 
average degrees of the layers are different (ranging from 3 
in the most sparse layer to 61 in the densest layer). In [3], 
the friendship ties of a group of students is mapped onto 
Facebook friends, picture friends and cohabitation, and 
these layers are shown to have different connectivity dis¬ 
tributions, hence different average degrees. In [33], the 
Indian airline and railway transportation networks are 
mapped onto two layers, representing two distinct modes 
of transportation between geographic locations. The de¬ 
gree distributions of these layers are then depicted, and 
it is observed that they have different degree distribu¬ 
tions (as well as different nearest-neighbor degree dis¬ 
tributions). In |34] . the international trade network is 
mapped onto 97 layers, each layer pertaining to one dis¬ 
tinct commodity. The connectivity patterns are differ¬ 
ent across layers. Although in this example the layers 
represent weighted networks, the assertion that the con¬ 
nectivity patterns are heterogeneous still holds. 

In the present paper, we consider heterogeneously- 
growing layers. First, we consider a simple two-layer 
system for expository purposes. We obtain n{k,i), the 
fraction of nodes with degree k in the first layer and de¬ 
gree £ in the second layer. We also use this result to 
find I{k), the expected layer-2 degree of a node whose 
degree in layer 1 is k. We solve the problem for the 
cases of preferential and uniform growth, separately. We 
demonstrate that the expression for i{k) is identical un¬ 
der these two settings. For the special case of homoge¬ 
neous growth, this result agrees with those of [29] . 

We then generalize both the preferential and the uni¬ 
form setups to M > 2 layers. Each incoming node es¬ 
tablishes pi, 132, ■ ■ ■, Pm links in layers 1,2, ...,M, re¬ 
spectively. For both uniform and preferential growth we 
obtain n(k), which is the fraction of nodes with vector 
degree k. That is, the fraction of nodes with degree ki in 
layer 1, degree /c 2 in layer 2, and so on. In all cases, we 
corroborate our theoretical findings with Monte Carlo 
simulations. 

In addition to heterogeneity of connectivity across lay¬ 
ers, this paper contributes to the literature by taking 
a step towards extending the analysis of the network 


growth process beyond the steady state, by considering 
arbitrary times. To that end, in Section [Vlj we focus on 
the uniform attachment model and analyze it in more 
detail. In particular, we study the temporal evolution 
of the joint degree distribution of a given arbitrary mul¬ 
tiplex network, with arbitrary number of layers, whose 
joint degree distribution is known. We find nt(k), which 
is the fraction of nodes with degree vector k at arbi¬ 
trary time t. Through several case studies, we verify the 
accuracy of our theoretical predictions through Monte 
Carlo simulations on multiple example topologies. We 
consider diverse topologies in order to examine how the 
initial conditions influence the evolution of degrees and 
to ascertain the accuracy of the predictions for different 
structures. Simulation results are consistently in good 
agreement with theoretical findings. 


II. MODEL 1: PREFERENTIAL 
ATTACHMENT IN TWO LAYERS 

The network is constructed by a set of nodes and two 
sets of links. There are two layers, each housing one 
set of links. This means that node x can have a set of 
neighbors in layer 1 and a different set of neighbors in 
layer 2. Similarly, the degree of node x in layer 1 can 
differ from its degree in layer 2. We denote the degree of 
node X in layer 1 by kx and in layer 2 hy £x- We denote 
by Nt{k,£) the number of nodes that have degree k in 
the layer 1 and degree £ in layer 2 at time t. The fraction 
of such nodes is denoted by nt{k,i). 

At the outset, the network comprises N{{)) nodes, with 
Li(0) links in layer 1 and ^ 2 ( 0 ) links in layer 2. The net¬ 
work grows by the successive addition of new nodes—one 
at each timestep. Each new node, upon birth, establishes 
Pi links to the existing nodes in layer 1 and P 2 links in 
layer 2. 

In this model, the probability that node x receives a 
link in layer I from an incoming node is proportional 
to kx- Similarly, in layer 2, the probability that node x 
receives a link from the newcomer is proportional to £x- 
At time t, when a new node is introduced the values of 
N{k,£) can consequently change. If a node with degree 
fc — 1 in layer 1 and degree £ in layer 2 receives a link 
in the first layer, its degree in layer 1 increments and 
turns into k, and N(k, £) increments consequently. If a 
node with degree k in layer 1 and degree ^ — 1 in layer 2 
receives a link in layer 2, its degree in layer 2 increments 
and consequently, N{k,£) increments. Moreover, if a 
node that has degree k in layer 1 and degree £ in layer 2 
receives a link in either of the layers, N(k,£) decrements 
consequently. Finally, each incoming node has degree 
Pi in layer I and degree P 2 in layer 2 upon birth, so 
each new incoming node increments A(/3I,/32) by one. 
The following rate equation captures the evolution of the 
expected value of Nt{k,£) by addressing the said events 
with their respective probabilities of occurrence: 
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Nt+i{k,e) = Nt{k,£) 

ik-l)Nt{k-l,e)-kNt{k,£) 


+ /3i 
+ P 2 


Li(0) + 2^it 

{£-l)Nt{k,£-l)-£Nt{k,£) 


Ski3i5e02- ( 1 ) 


^ 2 ( 0 ) + 2/32^ 

We can use this to write the rate equation for n{k,£). 
Using the substitution n{k,£) = {N{0) + t)n{k,£), we 
obtain 

[7V(0) + t] [nt+i{k,£) - nt{k,£)] + nt+i{k,£) = 
ik-l)Nt{k-l,£)-kNt{k,£) 


+ I3i 


+ (^2 


Li(0) + 2/3it 

i£-l)Ntik,£-l)-£Ntik,£) 


4/3i|5^/32- (2) 


U2(0) + 2/32t 

Now we focus on the steady state, that is, the limit 
as t —>■ 00 (the validity of this assumption is verified 
through simulations below). In this time regime, the 
time variations of n{k,£) vanish, and we also have the 
following simplifications 

N{0) + t _ 1 

~ “ 2 
1 

2 


lim /3i - , , 
t^oo Ui (0) -f ‘2f3\t 

i.„A. "<0) + ‘ 


( 3 ) 


t—>oo 1 ^ 2 ( 0 ) + 2/32^ 

So ([^ transforms into the following equation in the 


steady state: 


n{k, £) = 


{k — l)n{k — 1,1') — kn{k,£) 

2 

{£ — l)n(fc, £ — 1) — £n{k, £) 




k(3i ^102 5 


( 4 ) 

where we have dropped the t subscript in the steady 
state. This can be equivalently expressed as follows 


n{k,£) = 


k-l 
k + £ + 2 
24ft 4/3; 


n{k —!,£) + 


1 


k + £ + 2 


n{k, £ — i) 


' 2 + /3i+/32' 

This difference equation is solved in Appendix]^ The 
solution is 


n{k,i) = 


2/3i(/3i + l)fe(/32 + l) 




(2 + /3i+ft)fc(A: + l)t(t + l) 


k- Pi 


( 6 ) 


This agrees with the findings in [35], and those in |5^ 
[53] —in the special case of I3i = fi 2 - 

This result can be also simplified and equivalently ex¬ 
pressed in the following form: 

2(l + /3i+/32)! {k-l)\{£-l)\ 


n(k, £) = ■ 


{2 + k + £)l (/3i - 1)!(/32 - 1)! 

(fc-/3i+£-/32)! 
(fc-/3i)!(£-/32)!- 


( 7 ) 


Figure]^ is a depiction of the joint degree distribution 
for the symmetric case of Pi = £32 = 3. Figurej^pertains 
to the asymmetric case of Pi = 10 and P 2 = 10. Note 
that in both cases, the origin is at (/3i,/32), because the 
joint degree distribution is zero if for all k < Pi or £ < P 2 . 



FIG. 1. (Color online) The joint degree distribntion for pref¬ 
erential growth with di = 3 and P 2 = 3. Theoretical result 
is presented in §. Since the values decay fast in k and £, 
we have depicted the logarithm of the inverse of this func¬ 
tion, for a smoother output and better visibility. The joint 
distribution attains its maximum at fc = /3i and £ = P 2 - The 
contours are symmetric with respect to the bisector because 
Pi and P 2 are equal. Note that each axis begins at its corre¬ 
sponding value of P, so the origin is Pi, P 2 - 



FIG. 2. (Color online) logarithm of the inverse of the joint 
degree distribution for preferential growth with Pi = 10 and 
P 2 = 2. Theoretical solution is presented in The joint 
distribution attains its maximum at A: = /3i and £ = P 2 - The 
contours are skewed towards the k axis because Pi > P 2 - 
Note that each axis begins at its corresponding value of P, so 
the origin is Pi, P 2 - 

Let us verify ^ via simulations. For the first setup. 
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we consider /3i = 4 and /32 = 3, and plot n{k,£) for 
multiple example instances of k and £. The initial seed 
networks are identical star graphs with 6 nodes in both 
layers. The results are depicted in Figure]^ We observe 
that the simulation results reach the horizontal asymp¬ 
totes predicted by § relatively early on in the process. 
When the size of the network is iV = 300, the effects 
of the initial seed graph are already negligible. For the 
second example, we consider /3i = 2 and /32 = 3. For 
the initial seed network, we consider complete graphs 
with 4 nodes in both layers. The results are depicted in 
Figure]^ We again observe convergence to the theoret¬ 
ical predictions. Moreover, we observe that equilibrium 
emerges when the size of the network is 20 times the size 
of the initial seed network. Note that the asymptotes 
in Figure and those of Figure differ, because they 
pertain to different growth parameters. 



FIG. 3. (Color online) Comparison between the steady-state 
joint degree distribution and the theoretical prediction of § 
for growth under preferential attachment with growth param¬ 
eters /3i = 4 and P 2 = 3. The horizontal lines represent the 
steady-state asymptotes which accord with Q, and mark¬ 
ers represent simulation results. The initial seed networks 
are star graphs with 6 nodes in both layers. Several exam¬ 
ple value of {k,i) are chosen for illustrative purposes. The 
results are averaged over 100 Monte Carlo trials. 


We can see the entanglement between the two layers 
through the conditional average degree, which can be 
derived from (§. For the nodes who have degree k in 
layer 1, we can find their average degree in layer 2. Let 
us denote this quantity by £{k). In the calculations, the 
marginal degree distribution of layer 1 is invoked. This 
is obtained in [33 133]. and we also obtain it in Equa¬ 
tion (35) in Section IV —as a byproduct of our calcula¬ 
tions. We use this expression in the calculation of the 
conditional degree. To find the conditional degree, we 



FIG. 4. (Color online) Comparison between the steady-state 
joint degree distribution and the theoretical prediction of (|^ 
for growth under preferential attachment with growth param¬ 
eters /?i = 2 and P 2 = 3. The horizontal lines represent the 
steady-state asymptotes which accord with ([^, and markers 
represent simulation results. The initial seed networks are 
complete graphs with 4 nodes in both layers. Several exam¬ 
ple value of {k,l) are chosen for illustrative purposes. The 
results are averaged over 100 Monte Carlo trials. 


need to perform the following summation: 


m 


e e 


n{k, €) 
nk 


= E^ 


i 


/S, +/39+2\ 

2/3i(/3i-|-1)/32(/32-I-1) ( P^ + i ) ^fc-/3i-|-f-/32\ 

{2+/3i+P2)Hk+l)i{i+l) ('=+_(+=') V fe-/3i / 

k(k+l){k+2) 


/32(,d2 + l)(fc + 2) 

4- (2 + /3i+/32)(^+l) 


// 5 i +/ 32 + 2 \ /k — l3i-\-£—02\ 
\ /Si + 1 / V k—0i ) 


(k+£+2\ 
V /c+1 ) 



m 2 + 1 ) (^ 4 + 1 ^') 

(2 + /3i+/32) 


( 8 ) 


In Appendix [B| we perform this summation. The an¬ 
swer is 


m=^^^{k + 2). (9) 

In the special case of /3i = /32 = w, this reduces to 
, which is consistent with the previous result in 
the literature [29]. 

We now focus on the distribution of total degree (i.e., 
the sum of degrees in the two layers). Let us denote k+£ 
by q. The joint distribution of g, k is simply n{k, q — k). 
If we sum over all possible values of k, we get the dis¬ 
tribution of q. Note that k is at least /3i, because every 
incoming node has an initial degree of /3i in the hrst 
layer upon birth. Similarly, note that q — ki is at least 
j32- Taking these two into account for the summation 
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bounds, we have: 

2^1 (/?1 + l )/? 2(/?2 + 


n{q) =■ 


(2 + /3i + p2) 

9-/32 (q-Pi-P2\ 

V fc-/3i / 


X V 

fc(/c + l)(g - fc)(g - fc + 1) («+^) 
2/?i (A +1) A (A + 1) 

- n 


(2 + A + A)(9 + 2)(? + l)g(g-l) [izl 


We use the following identity: 

(%V) _ (/3i - 1)!(/?2 - l)!(g - 1) 


-/3i 

fe=/3i V/c-1/ 


(A+A-1)! 


( 10 ) 


( 11 ) 


This is proven in Appendix 0 Plugging this result 
into (ITol, we arrive at 


, X 2(a + A)(A + A +1) , o o\ l^<^\ 

"<’>= - ,(, + !)(,+ 2) - "(5-A-A). (12) 


This can be simplified and equivalently expressed as fol¬ 
lows 


(16) 


This difference equation is solved in Appendix The 
solution is 


n(fc, i') 




(1 + A + ■ 


(17) 


This result agrees with the findings in [35], and those 
in [28l |29]— in the special case of A = A- 

Figure represents the joint degree distribution for 
the case of /3i = A = 3. We plotted the logarithm of 
the inverse of this function for smoothness and visibility 
purposes. Figure [^depicts the joint distribution for the 
case of /?! = 10 and A = 2. Note that in both cases, 
since the joint degree distribution is nonzero only for 
k> j3i and £ > A, the origin is situated at /3i, A- 


where u{-) is the Heaviside step function. This is identi¬ 
cal to the expression for a single-layer network growing 
under growth parameter /3i -I- A- la other words, the ag¬ 
gregated network is tantamount to one which grows un¬ 
der the preferential attachment mechanism where each 
incoming node establishes /3i -I-A links to existing nodes, 
which is intuitively expected. 


III. MODEL 2: UNIFORM ATTACHMENT IN 
TWO LAYERS 


In this model, we assume that each incoming node 
establishes its link in both layers by selecting existing 
nodes uniformly at random. The rate equation ([^ in 
the case of uniform attachment transforms into 


iVt+i(fc,£) =iVi(fc,^)+A 


+ A ‘ 


A(0) -f t 


Nt{k - IJ) - Nt{kJ) 
N{Q)+t 

— + A/3iA,S2- (13) 


Using the substitution nk,e{t) = this becomes 


[N{0) -I- 1] [nt+i{k,£) - nt{k,£)] + nt+i{k,e) 


+ A 


+ A 


Nt{k-l,i)-Nt{kJ) 

N(0)+t 

Ntik,£-l)-Nt{k,£) 
'N(0) +t ■■ 


+ A/3iA/32, 


(14) 


which simplifies to the following difference equation in 
the limit as t —)■ oo, that is, the steady state: 

n p n{k-l,£)-n{k,£) n{k,£ - 1) - n{k,£) 

n{k, £) =A- - -h A-- 

+ A./3iA,,S2- (15) 



FIG. 5. (Color online) The joint degree distribution for uni¬ 
form growth with /3i = 3 and A = 3. Theoretical result is 
presented in fuf. Since the values decay fast in k and £, we 
have depicted the logarithm of the inverse of this function, 
for a smoother output and better visibility. The joint dis¬ 
tribution attains its maximum at k = Pi and ^ = A- The 
contours are symmetric with respect to the bisector because 
Pi and P 2 are equal. Note that each axis begins at its corre¬ 
sponding value of P, so the origin is di, A- 

Let us verify 0 via simulations. For the first setup, 
we consider di = 5 and A = 8. We plot n{k, £) for multi¬ 
ple example instances of k ad .^.The initial seed networks 
in both layers are rings with 10 nodes. The results are 
illustrated in Figure The simulation results visibly 
converge to the theoretical predictions. Moreover, equi¬ 
librium emerges relatively early, that is, even when the 
size of the network is 200, which is only 20 times the size 
of the initial seed network, the effects of the initial nodes 
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FIG. 6. (Color online) logarithm of the inverse of the joint 
degree distribution for uniform growth with /3i = 10 and 
/32 = 2. Theoretical solution is presented in Q. The joint 
distribution attains its maximum at k = and £ = j32- The 
contours are skewed towards the k axis because /3i > 02 - 
Note that each axis begins at its corresponding value of /3, so 
the origin is l3i,/32. 


are already negligible. For the second example, we con¬ 
sider the symmetric case f3i = P 2 = 10. For the initial 
seed network, we consider complete graphs with 15 nodes 
in both layers. The results are depicted in Figurej^ The 
results again verify the steady-state prediction, and con¬ 
firm that equilibrium emerges relatively early on in the 
process, namely, when the size of the network is 20 times 
the size of the seed network. 

To find the conditional average degree, we first need 
the degree distribution of single layers. This is found 
previously for example in [21[37]. The degree distribu¬ 
tion in the first layer is 


1 



Pi 

/3i + 1 


fc—/Si + l 


(18) 


We need to compute 

i i ^ 

ofe-/3i ot-P2 (k-px+l-P2\ 

Pi P 2 I fe-/3i ) 

^ (1 + /?1 + / 32 ) fe -/ 3 i +^-/32 + l 

/3i 

{Pi +P 2 + l)'=-/5i + l 4 - (1 + /?! + P 2 Y-f^^ 

(19) 

We have performed this summation in Appendix n 
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FIG. 7. (Color online) Comparison between the steady-s tate 
joint degree distribution and the theoretical prediction of (171 
for growth under uniform attachment with growth parame¬ 
ters /?! = 5 and (32 ~ 8. The horizontal lines represent the 
steady-state solution given by ( |17[ ), and markers represent 
simulation results. The initial seed networks are rings with 
10 nodes in both layers. The simulation results are averaged 
over 100 Monte Carlo trials. 


• (10,12) • (12,11) • (13,10) • (11,13) • (14,11) ■ (12,14) 



FIG. 8. (Color online) Comparison between the steady-state 
joint degree distribution and the theoretical prediction of 0 
for growth under uniform attachment with growth parame¬ 
ters Pi = P 2 = 10. The horizontal lines represent the steady- 
state asymptotes which accord with (171, and markers rep¬ 
resent simulation results. The initial seed networks are com¬ 
plete graphs with 15 nodes in both layers. Several example 
value of (fc, £) are chosen for illustrative purposes. The results 
are averaged over 100 Monte Carlo trials. 


The result is 

m = ^^ik + 2). ( 20 ) 

This is identical to ([^. For the special case of Pi = 
(32, this result (that the expected degree has the same 
expression under preferential and uniform attachment 
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schemes) is also obtained in [53]— see Equations (S.40) 
and (S.45) therein. 

Finally, to obtain the distribution of total degree, we 
undertake the steps similar to those in the previous sec¬ 
tion. We have: 


niq) = '^n{k,q- k) 

k 

ok-Pi oq-k-P2 (q-0i-l32\ 

E Pl P 2 I k-0i ) 

+ 


_ P 2 oh-Pio-k+Pi 7? - dl - /t 2 \ 

(1 + /3i-K/32)5-'9i-'9=+i \ k - P ^ ) 


{i+p, + -^<i^h-p2+i nj 

1 / P 1 + P 2 

Pi + P 2 \ ^ + Pi + P 2 ) 


q-Pl-P2 


( 21 ) 


This is similar to (18). This means that the degree dis¬ 
tribution of the aggregated network is identical to that 
of a uniformly growing network in which each newcomer 
forms Pi -\- p 2 links to existing nodes. 


IV. GENERALIZATION OF MODEL 1 TO M 
LAYERS 


Now let us consider the case of M layers, where 
the network possesses one set of nodes and M dis¬ 
tinct sets of links. Each incoming node establishes Pm 
links in layer m to existing nodes in that layer, where 
m G {1,2,. ..,M}. The initial number of links at layer 
m is denoted by Lm{0)- In this model, the degree of 
each node can be conveniently represented by a vector 
of length M. The degree vector of node x, denoted by 
ka;, stores the degree of node x in layer m as its m-th 
component. Let N (k) denote the number of nodes whose 
vector of degrees is k. Let n(k) denote the fraction of 
those nodes. If node x receives a link in layer m, then 
its degree will change to k -|- le™, where is the unit 
vector in m-th dimension, that is, it is a vector whose 
elements are all zero except its m-th element, which is 
unity. Let us denote the mth component of vector k by 
km- The rate equation for N{k) reads 

M 

Nt+i{k) =Nt{k)+ n 

m—1 

^ {km-l)Nt{k-em)-ikm)Nt{k) 

LmiO)+2Pmt 


Now we use Nt{k.) = and the following limit: 

N{0)+t 1 


lim /3„ 

t—^OO 


Lm(0) + 2Pmi 


(23) 


With these substitutions, we can rewrite (22) in the 
steady state as follows 


(km - l)n(k - em) - (km)n(k) 

M 

+ (24) 

m—1 

This can be rearranged and expressed equivalently as 
follows 



M 




km 1 




-n(k- e^) 


M 


2 + '^lm=l k'm m=l 
We first define the following auxiliary function: 




def (2 + X)m=l 


m = 


YC=l{km-l)\ 


l(k) 


(25) 


(26) 


Inserting this into (25) yields the following simplified 
equation 


M 


_ /c\ I 7 \| ^ 

P(k) = 5] (/.(k - e„) + 2 7 ,' n 


m—0 

M 


l\m=likm 1)! m=l 


= x : 0(k - e„)+2 n 


m =0 


Y\m=liPm 1 )! m= 


= 1 


(27) 


We define the M-dimensional Z-transform as follows: 

00 00 00 

^ — kl — k2—. .. — kM (28) 


ki—0k2—0 kM—0 


Plugging this into (27), we obtain 

M 

n 


$(z) = 2 


UZ=liPm-iy- ■ 

The inverse of this generating function is given by 


(29) 




UZ=iiPm - 1)! J Z=i 1 - E™=o (27rf) 


;\M ■ 


(30) 


(22) 



























1 


{2TTi)^ 

1 


nr?_ dz - _ - 

.j. , _ — 

= 1 ^ 2 ^m =0 


n 


(27rz) 


M 


M 

n 


^kr^-fSrr.- 


dZm ^ ^ 


n—0 


(27ri) 


M 


M 

n 

m—l 




= E E 

n=0 ri+...+rM=n 


dZn.J2 E 


^1 


V ri ,..., vm y 
n — Ori -\-...-\- rM—‘n ^ 


M 


n^i 


M 


E 

m—O 


, ■ ■ ■, ?'M/ \ki- Pi,..., km- PM 


( 31 ) 


Inserting this result into (30), we arrive at 


(2 + EEi/^™)! 


P(k) = 2 


nEi(^--i)! 


J2m=Akm - I3m) 

ki Pi,..., km Pm 


(32) 


Finally, from (|2^ we obtain 


i(k) = 


WZ^iikm-iy. (2 + EEi/3>n)! 
Y{l=i{Pm-l)\ (2 + EEiM! 

yM 


2 + 

' iL^m— 


'Em=likm- Pm) 
ki Pi,..., km Pm 


(33) 


Note that this expression holds for km > Pm, for all 
M layers. If any of the km values are smaller than the 
corresponding Pm value, then n(k) is zero. This is true 
because every node has an initial degree of Pm in layer 
m upon birth, and throughout the growth process, the 
degrees cannot decrease. 

Now let us consider some special cases in order to gain 
insight into the pattern (33) follows. For M = I, we 
have: 


„(t) = ltd)! E±4i; X 


k-p 

k-p 


(/?-!)! (2 + fc)! ■' 2 
It can be easily simplified and rewritten as follows: 
{k-l)\{P + 2)\ 2 2P{P + 1) 


(34) 


i{k) = 


iP - I)! {k + 2)\ P+ 2 k{k + l){k + 2) ■ 

(35) 


This agrees with the degree distribution of the prefer¬ 
ential attachment model for a single layer, obtained for 
example in [33 [53] . 

Now let us consider the special case of M = 2 and 
verify that (33) indeed reduces to ([^. For M = 2 we 
have: 


n{ki,k2) = 


(fci-l)!(fc2-I)! {2 + pi+p2)l 
(/?i - 1)!(/32 - 1)! (2 + fci+fc2)! 

2 (fci - Pi+k2- /32)! 
2 + Pi + P 2 {ki — Pi)\{k 2 — P 2 )' 


Note that the factor [2 + Pi-\- P 2 )\ that exists in the nu¬ 
merator of the second fraction on the right hand side 
simplifies into {1 + p + \ P 2 )\, due to the existence of 

the factor 2 + Pi + P 2 in the denominator of the third 
fraction. Applying this change, the result becomes iden¬ 
tical to the expression obtained in ([^. 

Finally, let us inspect the form of ( [Mj ) for M = 3. For 
three layers, we have: 

n{ki,k 2 ,k 3 ) 

{ki - iy.{k2 - l)!(fc3 - 1)! {2 + P 1 +P 2 + Ps)! 

(Pi - 1)!(P2 - 1)!(/?3 - 1)! (2 + ki+k 2 + ks)! 

2 (ki — Pi + k 2 — P 2 + ka — P 3 )! 


2 + Pi + P 2 + Ps (ki — Pi)!(k2 — /32)!(^3 — hY- 


(37) 


The pattern for general M is clear from these three cases. 


Let us verify the theoretical prediction in (33) via sim¬ 


ulations. We consider three layers, and an initial seed 
network with 5 nodes. In layer I, the nodes are situated 
on a ring. In layer 2, the topology is a star. In layer 3, 
we have a complete graph. The growth parameters are 
Pi = 2, P 2 = 3, and P 3 = 4. Figure [^depicts the results 
for multiple example (ki,k 2 ,k 3 ) triplets. In all cases, 
the simulation results visibly converge to the horizontal 
asymptotes predicted by (33). 


V. GENERALIZATION OF MODEL 2 TO M 
LAYERS 


Now we assume that there are M layers, all of them 
growing under uniform attachment. The rate equation 
reads 


M 


Nt+i{k) =Ntik)+ n 


m—l 


M 


+ E 


771=1 


Nt(k-em) - Nt{k) 
Ni0)+t 


X 


(36) 


(38) 
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Let us define 


B=l + i:l=rlim 

(^f Pm 


(41) 


Taking the generating function of two sides of (40), we 
get 




^ M M 

g n +^(^) 

m—l m—1 


(42) 


FIG. 9. (Color online) Comparison between the steady-state 
three-layer joint degree distribution and the theoretical result 
of (331 for the case of preferential attachment. The horizontal 
lines represent the steady-state asymptotes predicted by (331, 
and markers represent simulation results, which are averaged 
over 100 Monte Carlo trials. The initial seed networks have 
5 nodes, and the topologies in the three layers are ring, star, 
and complete, respectively. The growth parameters are /3i = 
2, dz = 3, and da = 4. Several example values of (fci, k 2 , fcs) 
are chosen for illustrative purposes. 


This can be equivalently expressed as follows 


V'(z) 


1 

B 


M 

n 


M 


1 - 

m—l 


(43) 


For rii(k), this transforms into the following recurrence 
relation in the steady state 


M 


M 


n(k) = /3n 


nt(k - e^) - nt(k) 


(39) 


This can be rearranged and recast as 


M 


n(k) = 


n 




M 




n((k- e^) 


m—l 




(40) 


The inverse of this generating function yields the desired 
degree distribution: 


n(k) 


1 1 

B {2TTi)M 


TT 

J_ m m 
m—l 
M 

^ ^ Qm^rn 
m—l 


(44) 


We can perform the integration by undertaking the 
following successive steps: 
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(27rf 


M 


M 

n 


-/3„ fc™-l 

m 


dZri 


M 


i-E 


Qm ^ri 


(27rf 


M 


M 


oo r “in 

-1 




771=1 


71=0 


771=1 


(27r?) 


M 


M 

n 

771=1 




= E E 


n=0 ri + ...+rM—IT' ’ * ’ ’ ’g— 


^dZm Y. E 




EEl(fcm -/^™) 

^1 — A, • ■ ■, kM — Pm 


M 

n 


^km — ^ 


M 


n=0ri + ...+rM=n ^ ^ s=l 

M 




(45) 


Plugging this result into (44), we get 


n(k) = ^( ^ 


n 


.krr^-^ 


‘. (46) 


Using the definition of B as given in (411, this can be 
expressed equivalently as follows: 


M 


( 'I TT «/c„-/3„ 

Ifei —/9i,...,feM —/SmI 11”" 

«(k) = -- w (47) 


(i + EEi/3-) 


l + E“=l(fem-/3r„) 


Note that, similar to the case of preferential growth, this 
expression holds for km > Pm, for all M layers. If any 
of the km values are smaller than the corresponding Pm 
value, then n(k) is zero. 

Now we consider the special cases of M = 1,2,3 and 
simplify this expression to gain more insight into the 
pattern it follows for general M, for convenience of in¬ 
terpretation. 

For M = 1, we have 


n(A:) = - 


l + (fc-/3) ' 


(1 + / 3 ) 


(48) 


The binomial coefficient equals unity, and the result be¬ 
come identical to Equation (18), which is also obtained 
for example in 


Now we focus on the special case of M = 2, and con¬ 
firm that it agrees with 0 - Note that the multinomial 
coefficient in (47) becomes the ordinary binomial coeffi¬ 


cient in this case. So for M = 2 we have: 


n{ki,k2) = 


(ki-pi+k'2-P2\ oki-Pi qk2/32 

I fcl-/3l ) Pi P2 

{l + Pi+ ’ 


(49) 


r 


which is identical to Equation (18). 


Finally, let us consider the case of M = 3. We expand 
the multinomial coefficient into factorials to render the 
pattern for general M more apparent. We have: 


(ki — dl + ^2 — fe + ^3 ~ fe)! qki-l3i 


n(ki,k2,k3) = 


- Pl)Kk2 - p2)Kk3 - PsY- 


“'’"da 


-A 


(1 + di + da + da) 


1 -\- ki — 01 k2 — 132 - 


-03 ■ 

(50) 


Let us verify the theoretical result of (47) via simula¬ 


tions. For the initial seed graph, consider 8 nodes which 
form complete graphs in layers one and three, and let 
them form a star graph in layer 2. Let the growth param¬ 
eters be Pi = 3, P 2 = 4, and P 3 = 2 . Figure [T 0 | depicts 
the results for multiple example (ki,k 2 ,kz) triplets. In 
all cases, it is visible that the simulation results converge 
to the horizontal asymptotes predicted by (ItT 


VI. 


TEMPORAL EVOLUTION OF THE JOINT 
DEGREE DISTRIBUTION 


We now move beyond the steady-state analysis and fo¬ 
cus on the temporal evolution of the structure of multi¬ 
plex networks. Consider a given multiplex network with 
V(0) nodes, whose joint degree distribtuion is given by 
no(k). How would n(k) evolve if new nodes are added to 
the network? In this section, we seek a time-dependent 
solution for nt(k). We only consider the case where at¬ 
tachments are uniformly at random in each layer. Let 
us focus on the evolution of iVt(k). At time t, upon the 
introduction of a single new node, we can rewrite the 
time-continuous analog of Equation (22) as follows: 
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FIG. 10. (Color online) Comparison between the steady-state 
three-layer joint degree distribution and the theoretical result 
of (|47|) for the case of uniform attachment. The horizontal 

(|47l), 


lines represent the steady-state asymptotes predicted by i 
and markers represent simulation results, which are averaged 
over 100 Monte Carlo trials. The initial seed networks have 
8 nodes, and the topologies in the three layers are complete, 
star, and complete, respectively. The growth parameters are 
/?! = 3, /32 = 4, and fis = 2. Several example values of 
(fci,fc 2 ,fc 3 ) are chosen for illustrative purposes. 


51Vt(k) 

dt 


M 


m—l 

M 


+ ^ /3m 

m—l 


Nt{k-e^)-Nt{k) 

N{0)+t 


(51) 


Let us define the following generating function: 

vl/,(zi, ..., ^m) =' E ■ • ■ E (52) 

All 


We now multiply both sides of (51) by Zi 
and sum over all kj. The result is 




5T't(zi,...,ZM) 

dt 



2=1 


M 




- ^ -(z-i 

N{Q)+e * 


l)T't(zi,... ,zm). 


(53) 


We can rearrange the terms and express this equation 
equivalently as follows: 

..., zm) ^ ~ l)lE't(2i) ■ ■ ■ 5 zm) _ TT 

di ^ iVYOf+t “IP* • 

' ' i=l 

(54) 


This is a linear first order equation |38j which can be 
solved by the multiplication of both sides by the follow¬ 
ing integration factor: 


» M 

fj-t{zi,Z m) = exp / - E 


- 1 ) 


dt 


= exp - ln(Ar(0) + ^ - 1) 

= (iV(0)(55) 
Using the integration factor, the solution is given by 


M 


'^t{zi, - ■ ■ iZm) 


-1 


= k' 


ft M 

z~^' fj,dt + C 

i=l 


= (/V(0)+t) 


1-EE/3*(P'-1) 


c 


i-EUA(P^-i) 


C{N{0)+t) 


Jl—1 

E"i/3.cr'-i) 


(56) 


In this equation, C is a constant in time, but can depend 
on zi,..., Zm- It must be determined in order to satisfy 
the initial conditions. Note that the generating function 
of the joint degree distribution of the initial network is 
known. So we can use tE'o( 2 i,..., zm) in order to identify 
C. Setting t = 0 in (56 1 , and rearranging the terms, we 
arrive at 


C =^o{zi, ...,zm)x /V(0)- ^*=1 ftp -1) 

.iv(0)i-^"iftp'-i). 

(57) 






Finally, plugging this into (56) and rearranging the 
terms, we obtain the generating function of the joint 
degree distribution at arbitrary times: 


4't(zi,...,ZM) = {N{0)+t) 

+ 410(2:1,..., Zm ) X 1 -|- 

t 


HZ 


M -Pi 


i=l 


i-EEaP-pi) 




W(0), 


M 


n ivi 


- 1 ) 


(58) 


Now we need to take the inverse of this trans¬ 
form. First let us take the inverse transform of 

denote 1 -f Em=i /3m 

by B, and we denote ^ by gm, for brevity of notation. 
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Using straightforward properties of geometric series and 
multinational expansions, we have: 




B \ — . QjZ- B 

rn^l i 


m ri-\-...-\-rM—'f^ ^ ^ s—1 

The inverse transform of the expressions is given 
by the delta function 6 ^ . Using this result, we get 


M 


E 


m 


^ "V r, + .^rM=ra V^l> • ■ ■ > ru) 


M 


(60) 


Applying the delta functions, all the sums reduce to a 
single nonvanishing term: 


°'“"bU. 


M 


The multiplication by Jl^i U merely results in 
shifts, that is, every kj. changes to — Pr- So we obtain 


M -Pi 


n ivi 

i = l^: 


- 1 ) 


U E™rifc™-A„ 


B\ki-Pu---, kM - PmJ U 5'=^-'’" 


(62) 


Now let us take the inverse transform of 

E at /q ^ — 1 
i = l Pi^i 

G '2 =‘ n + 1 • Using elementary proper¬ 


ties of the exponential function, we have 

M , 

G.=n(i 


M 


= n 


m—l 


t 

W) 


In ( 1 + 


iV(0) 


M 00 

HE- 


In ^1 + jv(o) 


I 


(63) 


m—l r—1 


The inverse transform of the terms yield delta func¬ 
tions, which helps us eliminate the sums. The result is 


M 


Go 


lift- 


In (^1 + Ar(o)) 


n k„ 


k I 


(64) 


Plugging the expressions for the inverse transforms 
found in (62) and (64) into (58), and noting that the 


multiplication in the z domain is equivalent to convolu¬ 
tion in the k domain, we arrive at the following expres¬ 
sion for the inverse transform of , Zm), which 

is equal to Nt{zi,..., zm) by definition. Dividing the re¬ 
sult by the number of nodes at time t, which is N{0) + t, 
we obtain the joint degree distribution as: 


1 / k - B \ 

n 1 \ ft I Z^m=l 'ft”! P”! 1 TT Pr 

= B (*, - A, ..., - Pm) n Bftftft 


+ 1 


-B 


M 


N{0) 


1 / t 
B\^ N{0) 


U no(Ci,. ■ ■ ,Cm) P 

Cl.■■■,5 m ”1=1 


^ftft + N{0) 


km-ir. 


-B M 


n bp- y. ■■■ y. 


{kjYi Cm)- 


Cm Pn 


m=l 4i=/3i C«=/3 m PmJ Bir^{km Cm)! 


M 

n 


[in + N{0) 


Okm-ir, 


(65) 


r 


The first term on the right hand side of (65) is inde¬ 


pendent of time, and is the steady-state solution which 
was obtained in (47). The second term characterizes 


the effect of the initial network. Note that the factor 

/ \ —B —1 

M -|- ) makes this term vanish in the long-time 

limit. This is expected, because in the limit as t —?> oo, 
almost every node is among those who were subsequently 
appended to the network, not those who belonged to the 
initial network. The last term on the right hand side is 


also transient and vanishes as time goes to infinity. At 
t = 0 , note that the first and third terms on the right 
hand side of (65) cancel out, and only the second term 


remains. In the second term, all the individual terms 
in the sum are zero (because the logarithm in the sum¬ 
mand become log[l-|- 0 ], which is zero), except one term, 
which pertains to Cm = kmSm (and for the logarithm, 
0 ° = 1 ), which means that the result correctly reduces 
to no{ki,...,kM)- 

For the special case of M = 2, we get the following 
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expression for the joint degree distribution of two layers: 


nt{ki,k2) = 


fci - + fc2 - , 

fci - j3i 




(/3l + /32 + l)l+fci-/3i+fe2-/32 


+ 1 


^ ^ -(/3i+/32 + 1) 


iV(0) 




-5l fik2-^2 


In (^1 + Ar(o)) 


^ 1 —^ 1+^2 —^2 


^15^2 


- 1 + 




NiO) 


ki k2 

(Pi+ I32 + 1V-P^-^^ ^ ^ 


(fci-ei)!(fc 2 - 6 )! 

- ,dl + ^2 - /32 

6-/3i 


In 1 + 


N{0) 


ki — ^i-\-k2 — ^2 


iPi + P 2 + l)fi+42(fci — ^i)!(fc2 — C 2 )! 


( 66 ) 


In the limit as t —t 00 , the second and third terms on the right hand side vanish, and the first term prevails, which 
coincides with the result obtained in 


Now let us verify the theoretical prediction for the 
time-dependent joint distribution via simulations. The 
first setup we consider is as follows. Consider 50 nodes 
with identical connectivity in both layers at time t = 0, 
forming a ring. For the growth mechanism, consider 
,01 = 3 and /32 = 2. Figure [IT] depicts nt{ki,k 2 ) for some 
example test values for ki and fc 2 - It can be observed 
that simulation results are in good agreement with the¬ 
oretical predictions. 



ate neighbors to the left and b to the right. Then, each 
nonexistent link is established independently with prob¬ 
ability p. We denote such a network by SW{N{0),b,p), 
where SW denotes small-world. For the second simula¬ 
tion setup, we consider a SkF(200, 0.05, 2) for the first 
layer and a S'kF(200,0.01, 2) for the second layer. We 
set the growth parameters to be /3i = 3 and /32 = 4. The 
results for several test cases of fci and ^2 are depicted in 
Figure Good agreement is observed between simula¬ 
tion results and theoretical predictions. 



FIG. 11. (Color online) Temporal evolution of the joint de¬ 
gree distribution for two layers. The growth parameters are 
01 = 3 and 02 = 2. Th e solid curves represent theoretical 
prediction given by ( 661 . The markers represent simulation 
results, averaged over 100 Monte Carlo trials. 


For the second example, we consider the initial net¬ 
works to be small-world graphs [31]. We construct the 
initial networks in the following way. Consider a network 
of N{0) nodes. Suppose that the nodes are situated on 
a circle and then each node is connected to b immedi- 


FIG. 12. (Color online) Temporal evolution of the joint de¬ 
gree distribution for two layers. The initial networks are 
S'1F(200,0.05,2) and SVF(200, 0.01, 2), respectively, as dis¬ 
cussed in the text. The growth parameters are 0i = 3 and 
02 = 4 . Solid lines represent theoretical predictions, given 
in ( |66[ |. Markers represent simulation results, averaged over 
100 Monte Carlo trials. 

We can also illustrate the evolution of the joint distri¬ 
bution by depicting it at different instants of time. To 
that end, since at most three dimensions can be plotted 
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and nt{ki,k 2 ) has four dimensions, we need to discard 
one dimension to enable us plot the distribution at dif¬ 
ferent time points. To that end, we need to fix one di¬ 
mension. Since the index of layers are arbitrary, without 
loss of generality, we fix ki. For the simulation, consider 
the following setting. The initial network in layer 1 is a 
^^^(SOO, 2,0.02) and that of layer 2 is S'1F(300,4, 0.04). 
The initial number of connections of incoming nodes are 
(3i = A and (^2 = 2. We plot nt{ki, ^ 2 ) at different time 
steps as a function of fc 2 , having fixed ki at the example 
value of ki = 8. In other words, we plot nt(8,/c 2 ). The 
results are depicted in Figure [T^ 




FIG. 13. (Color online) The joint degree distribution 
nt{S,k 2 ) as a function of for different times. The ini¬ 
tial number of connections of incoming nodes are /li = 4 and 
/32 = 2. The initial network in layer 1 is a 5'IF(300, 2, 0.02) 
and that of layer 2 is SW (300, 4,0.04). The shades represent 
theoretical predictions, given in ( 661 . Markers represent sim¬ 
ulation results, averaged over 100 Monte Carlo trials. Shades 
are chosen instead of linear plot for better visibility. 


As another example, consider a ring with 100 nodes 
to be the initial network in the first layer, and let the 
nodes form an Erdos-Renyi (henceforth ER) network [30] 
in which the probability of existence of each link is 0.03. 
The results for the joint degree distribution with ki fixed 
at 2 is depicted in Eigure[T4| In other words, Eigure[T4| 
depicts nt( 2 ,k 2 ) as a function of ^2 at multiple instants 
of time. As time progresses, a peak at ^2 = h emerges. 
This is because as more new nodes enter the network, 
the number of newcomers who have degree (32 increases 
and they constitute the majority of the network at long 
times. This causes the initial lump in Figure [14] —which 
pertained to the initial network—to lose its mass and 
become smaller. The mass moves towards k 2 = /32 and 
construct a peak there. 

For the next example, we consider two layers with 
Barabasi-Albert (henceforth BA) topology |3Tj . The ini¬ 
tial networks have 300 nodes, and the initial connec¬ 
tivity of incoming nodes are 2 for both layers. So in 
the resultant network, the majority of nodes have de¬ 
gree 2, and the degree distribution has a power-law tail 
with exponent 3, as is expected in BA graphs. After the 
BA network with 300 nodes has been constructed under 


FIG. 14. (Color online) The joint degree distribution 
nt(2, ^ 2 ) as a function of ^2 at different instants of time (i.e., 
we have fixed fci = 2 to enable graphical illustration). The 
growth parameters are /3i = 2 and /32 = 8. The initial net¬ 
work in the hrst layer is a ring with 100 nodes, and in the 
second layer the topology is an ER network with the link 
existence probability equal to 0.03. The solid lines represent 
theoretical predictions, given in ( 661 . Markers represent sim¬ 
ulation results, averaged over 100 Monte Carlo trials. Note 
that as time progresses, a peak at k 2 = /32 emerges, which is 
intuitively expected. 


the said mechanism, we start the growth process with 
growth parameters /3i = 1 and (32 = 25. This means 
that the second layer will be more densely connected as 
compared to layer 1. We fix the layer-1 degree at fci = 3 
and plot nt(3, ^ 2 ) as a function of at different instants 
of time. The results are presented in Figure Note 
that the initial distribution has a peak at ^2 = 2, which 
is expected, because as mentioned above, the majority 
of nodes in the initial BA graph have degree 2 in both 
layers. As the growth process proceeds, the center of the 
peak moves towards higher values of k 2 - This is because 
the newcomers have initial degree ^2 = 25 in the second 
layer, and as time progresses, the number of such nodes 
increase. 

Let us also investigate the effect of sudden change in 
the growth parameters through one example. Consider 
a random recursive tree [4^ , which is constructed as fol¬ 
lows. We begin by a single node, and then new nodes 
enter the system sequentially, and each newcomer ran¬ 
domly chooses one existing node and connects to it. The 
result is a tree, called a random recursive tree (here¬ 
inafter RRT). We choose the initial networks in both 
layers to be RRTs with 100 nodes. We then apply to this 
substrate a growth process with /3i = 1 (which would be 
a continuation of the growth mechanism that gave rise 
to the RRT in the first layer) and (32 = 10. This means 
that in the second layer, the growth mechanism is un¬ 
dergoing a regime shift, as the number of connections 
per incoming node abruptly jumps from 1 to 10. The 
results for the temporal evolution of 714 ( 1 ,^ 2 ) is plotted 
in Figure We observe that a new peak emerges at 
k 2 = 10, and the initial lump that peaked at ^2 = I 
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FIG. 15. (Color online) The joint degree distribution 
nt(3, ^ 2 ) as a function of k 2 at different instants of time (i.e., 
we have fixed fei = 3 to enable graphical illustration). The 
growth parameters are /3i = 1 and /32 = 25. The solid lines 
represent theoretical predictions, given in ( 661 . Markers rep¬ 
resent simulation results, averaged over 100 Monte Carlo tri¬ 
als. The initial networks are BA networks in both layers, as 
described in the text. 


the third layer. For the growth mechanism, we consider 
/3i = I, (32 = 4, and = 5. The results for several 
example triplets of (fci, ^2, fca) are depicted in Figure [17 


. simulation results are visibly in good agreement with 
theoretical predictions. 



at the outset gradually loses its mass to the new lump 
concentrated at ^2 = 10, which embodies the nodes that 
are introduced to the network in the second phase of the 
growth process. 



FIG. 16. (Color online) The joint degree distribution 
nt(l, ^ 2 ) as a function of /c 2 at different instants of time (i.e., 
we have fixed fei = 1 to enable graphical illustration). The 
growth parameters are /3i = 1 and p 2 = 10. The solid lines 
represent theoretical predictions, given in ( | 66 [ ). Markers rep¬ 
resent simulation results, averaged over 100 Monte Carlo tri¬ 
als. The initial networks are RRT networks in both layers, 
as described in the text. 

Let us also consider one example to verify accuracy 
of the theoretical prediction of the joint degree distri¬ 
bution for more than two layers. Let us consider the 
following setting for the initial network: a BA network 
of 100 nodes with parameter 1 (that is, the initial num¬ 
ber of links each incoming node establishes upon birth is 
1) in the first layer, a BA network with parameter 2 in 
the second layer, and a BA network with parameter 3 in 


FIG. 17. (Color online) The joint degree distribution 
nt{ki,k 2 ,k 3 ) for several example triplets {ki,k 2 ,k 3 ). The 
growth parameters are /3i = 1, /I 2 = 4, and da = 5. The 
solid lines represent theoretical predictions, and the markers 
represent simulation results averaged over 1000 Monte Carlo 
trials. The initial networks are BA networks in all three lay¬ 
ers, as discussed in the text. 


VII. SUMMARY AND OPEN PROBLEMS 

In this paper we focused on heterogeneous growth of 
multiplex networks. We considered both the preferen¬ 
tial and uniform growth mechanisms. We analyzed the 
problem for M > 2 layers and obtained the joint degree 
distributions. For the case of uniform growth, we ana¬ 
lyzed the problem for arbitrary times, and obtained the 
joint degree distribution for arbitrary number of layers 
at arbitrary times, for given initial conditions. We veri¬ 
fied our findings through Monte Carlo simulations, and 
observed that the theoretical predictions are remarkably 
accurate. 

An immediate generalization of this work is to con¬ 
sider nonzero couplings in the preferential growth mech¬ 
anism, so that the link reception probability of a node in 
each layer depends on its degrees in all layers, and then 
find the joint degree distribution via the rate equation 
approach (nonzero couplings are envisaged in p5H5D] . 
but the rate equation approach remains unsolved and 
no closed-form solution for the joint degree distribution 
exists in the literature). For example, in the two-layer 
case, the connection kernel would depend on the degrees 
of existing nodes in both layers. Consider an existing 
node X, who has degree kx in layer 1 and degree £x in 
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layer 2 at time t. Then its probability of receiving a link 
in layer 1 from the incoming node would be proportional 
to giikj. + gi 2 (-x I and its probability of receiving a link 
in layer 2 would be proportional to g 2 ikx + <722-^2:• Our 
solution presented in this paper pertains to the special 
case of gii = (722 = 1 and gi 2 = 521 = 0. The rate 
equation corresponding to the general case (counterpart 
of Equation Q, with all four gijS nonzero, would be as 
follows: 


Nk,i{t + 1 ) = Nkj{t) + Sk0iSgp^ 


The following holds 

K + t + 2 n 


^ n{k,i-l)=& -- ^^{k,e.)[k + ^ + 2)\mk,i-i. 


k + (. + 2 

Plugging these into , we can recast it as 

{Pi + /32 + 1)! 


Tfiki = ruk-if + mkA-i + 2 


(/?! - IW 2 - 1)! 


(A3) 

5k,Pi5l^P2- 

(A4) 


■Pi 


■ P 2 


-Pi 


gii{k - l)+gi 2 ( 


Nt{k-1,£) 


-P 2 


giiLi 

( 0 ) + 312(62(6) -I- 

321 A; -1-522(6- 1) 

2 .{giiPi + 312^2)^ 

Nt{k,£-l) 

52 iLi(0 ) 4 

322-62(0) -I- 

3116: + 3126 

2 {g 2 iPi + 322/32)^ 

Nk,e{t) 

5 iiLi( 0 ) 4 

312-62(0) -I- 

321^-1-3226 

2 {giiPi + gi2P2)t 

Nk,e{t) 


g2lLl{0) + 322 -^ 2 ( 0 ) + 2{g2lPl + g 22 P 2 )t 


(67) 


It is straightforward to simplify this equation in the 
steady state, but solving the resulting difference equa¬ 
tion is more intricate than the one studied in this paper. 

Furthermore, a related quantity which is ubiquitous in 
the studying of epidemics and various diffusion processes 
over networks is the nearest-neighbor degree distribu¬ 
tion [33H1S]- Let us denote this quantity hy pi{k,£\q,r). 
For a node with degree k in layer 1 and degree £ in layer 2, 
if we select one of its neighbors in layer 1 randomly, then 
Pi{k,£\q,r) would be the probability that this neighbor 
has degree q in layer 1 and degree r in layer 2. Similarly, 
P 2 {k, £\q, r) would be the probability that a randomly se¬ 
lected layer-2 neighbor of a node with degrees k,£ will 
have degrees <7, r. This quantity is essential for studying 
dynamics on networks, and to obtain it for multiplex 
networks, one also needs exact expressions for the de¬ 
gree distributions—which the present paper focused on. 


Now define the Z-transform of sequence rrik^e as fol¬ 
lows: 




k i 


rrik/ = j i’{z,y)z'" ^y^ ^dzdy. 


(A5) 


Taking the Z transform of every term in (A4), we ar¬ 
rive at 


il){z,y)=z ^'ti){z,y)+y V(^,2/) 

2 (/^ i +/^2 + l )! -01 -t 

{Pi-mP2-l)\ ^ 


(A6) 


This can be rearranged and rewritten as follows 


ip{z,y) 


2 


1 — z ^ — y ^ 


(/?l+^2 + l)! -01-02 

{Pi - l)\{p2 - 1)! ^ 

(A7) 


The inverse transform is given by 

2{Pi+P 2 + 1)\ ff z^~^^~^y^~^^~^dzdy 
(/3i - 1)!(/32 - 1)! / / (-47r2)(l - z-i - y-i) 
2{Pi+P2 + l)\ ff 
{Pi - 1)\{P2 - 1)! / f {-^TT'^){zy - z-y) 

2{Pi +P 2 + 1)! f f z^-f^^y^-^^dzdy 
{Pi - 1)\{P2 - 1)! / f (-47r2)(y -l)[z- 

(AS) 

First we integrate over z. We get 


Appendix A: Solving Difference Equation ([^ 


1 


We need to solve 
fc — 1 

^ ^ k + £ + 2 ^ ^ k + £ + 2 

2Sk/3i Sei3^ 

2 + Pi + P 2 

We define the new sequence 

def {k -£ £ -£ 2)\ 


i{k, £ — 1) 


mkg = 


(A:-!)!(£- 1)! 


n{k, £). 


(Al) 


(A2) 


_ 2{Pi + P 2 -\-xy, f (y-i) y ^y 
{Pi - 1)1{P2 - 1)1 J (27ri)(y-l) 
2{Pi+P 2 + 1)\ f y^-f>^+^-P^dy 
~ {pi-l)\{P2-iy.J (27ri)(2/-l)fc-/5i+i- 

(A9) 


Now note that the residue of 


f{y) f 

7-— tor positive 

(j/-l)” 


integer equals ^-7,-, where the numerator denotes 

(n- 1)! 
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the (n — l)th derivative of the function /(y), evaluated 
at y = 1. Also, note that the m-th derivative of the 

m\ 


function y", for integer n and m, equals 
Combining these two facts, we obtain 


(n — to)! 


y 


f/^l-+02+2\ fk-Pi+l-P2\ 

m = ^ ’ (Bi) 




2(/3i+^2 + 1)! fk — Pi+£ — l32\ 
Using (|A2|), we arrive at 


n{k,£) = 


2(^i+/32 + 1)! 1 (" 


(/li - 1)!(/32 - 1)! + mi + 1) {'^+1+^) 

(All) 

This can be equivalently expressed as follows: 

2dl(ft + l)fe(d2 + 1) A - ft + £ - fe 


n{k,l) = 


(/3i+/?2 + 2)fc(A: + l)df + l) ('=+|f) V k-Pi 


(2 + /?! + /?2) 

We use the following identity: 


74y = (n + l) /\-(l-t)—dt, (B2) 

U) JO 

to rewrite the binomial reciprocal of the coefficient as 
follows 

^^ = (fc + £ + 3) (B3) 

\ t ) Jo 

Also, from Taylor expansion, it is elementary to show 
that 


(A12) 


Si{x,n) ^ 


def \ ^ 'KK. / ^ 
X 

n 


Appendix B: Performing the Summation in ^ 
We need to calculate 


(l-x)”+i' 

This identity will be used in the steps below. 


(B4) 


Plugging (B3| into (B5), we have 
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m = E 


/32(/32 + 1) 


(2 + + /32) 




('■ A ? r 0 '0 i‘ ‘'<‘ - 


(2 + /?! + /I 2 ) 

/?2(/32 + 1) //3l + /32 + 2 


(2 + /3i + /I 2 ) V /3i + 1 

^2{/32 + 1) //3i + /32 + 2 

(2 +/3i +/32) V /3i + 1 / Jo 

/32(/32 + 1) //3i +/32 + 2^ 


^ ^ (1 - E(^ + ^ + 3)t''+^+^ 


dt 


(2 +/3i +/32) V /3i + 1 / Jo 

m /32(/32 + 1) (Pi + ^2 + 2^ '■i 


d 
dt 


E*" 


fc+l+3 (k - Pi + I - p2 

k- Pi 


dt 


(l-ty+H 


k-\-2±—k—2 


^3+/3i+/32 


k-i3i+e-P2 fk - Pi + ^ - P 2 

k-Pi 


(2 + Pi-\-P2) V /3i + 1 
P 2 {P 2 + 1) f Pi + P 2 + “2: 


(1 - ty+H-'^-^ 

^ ' dt 


^3+/3i+/32 


^fc-/3i 


(l_t)fe-/3i + l_ 


dt 


{2 + Pi + P 2 ) V /3i + 1 / Jo 

/?2(/32 + 1) f Pi + P 2 +2 


\l_t)fe+2i-fc-2^ 


^fc+/32+3 


(2 + /?! + /I 2 ) V /?! + 1 

_ ,52{/32 + 1) //3i +/32 + 2\ 

(2 + /?! + P 2 ) \ /3i + 1 / 

P2{P2 + l) /^/^i+^2 + 2\ r ^i!/32! . . /3i!(/32 + 1)! 

- (2 + A + A)l A + i j + * + 3) (ft + + 1), - <1 + '*■ + *> (A + & + 2)!j 




(1 - ty-di+i 

^ ^ (1 - [fc + /?2 + 3 - (1 + /3i + /32)i] dt 

{k + P2 + i) f {1 - y'^^y^dt - {1 + Pi + P 2 ) f {1 - t)^^y^~^^dt 

JO Jo 


P2{P2 + P)PPP2' 


(Pi + P 2 +2 

{2 + Pi+ /^2)(1 + Pi+ P2y- V /3i + 1 

/?2 


[{k + /32 + 3) — {P 2 + 1)] 


Pi + 1 


ik + 2) 


(B5) 


Appendix C: Proving the Identity Given in (111 


We u se the property of the Gamma integral given 
in (B21 in order to deal with the binomial coefficient 


that is in the denominator of the summand. We have: 


E 


V k-Pi ) 


k \k-l) 

_ (q — Pi — P 2 


k-Pi 


= (g - 1 ) (1 + 

= ((7-1) - ty^-^ 

JO 


dt 


q-di-02 


dt 


jB2l 


(<7-l) 


(/32+/3i-1)(%+/_E) 

(/?i-1)!(/32-1)!((7-1) 

(/32 + /?i-1)! 


(Cl) 
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Appendix D: Solving Difference Equation (16l After inserting the expressions for qi,q 2 from (D2), this 

becomes 


Let us repeat the equation we need to solve for easy 
reference 




n{k,£) = 


qk-Pl (k-Pi+e-P2\ 

P2 I k-l3i ) 


(Dl) 


Let us define the following quantities from brevity: 


(1 + /3l + /32)'=-/5 i+^-/52 + 1 • 

Appendix E: Performing the Summation in 

We need to perform the following summation: 

^2 V fc-/3i 


(D6) 


m = 


+ l)fe-/3i + l ^ (k-l3i+t-l32\ 




def /3l 

“ I + P 1 +P 2 

def _ h 


(/3i + /32 + l)'=-/5i+i ^ (1 + /3i + 

(El) 


(D2) 


Let us denote k — j3i by k' and £ — /32 by Also 
let us denote by a;. De need to evaluate the 

following sum: Let us use (B4| 

Taking the Z transform from both sides of (Dl), we define 


get 


'fp(z,y) = qiz ^'fpiz,y)+ q 2 y V(^, 2 /) + 


5'i(a;,n) ^ 


Z-Ply-P2 


m 


def \ ^ ^ ^ \ 

n J (1 — a;)"+^ 


(E2) 


1 + /?i + /32 ' We have: 
(D3) 


This can be rearranged and recast as 
1 


k' + e 


ip(z,y) = 


^-Piy-P2 


= /32X~^ Si{x, k') + X £'x^ 


1 - qiz~^ - q2y~^ 1 + /3i + /32' 

This can be inverted through the following steps 
1 


(D4) 


_i I k'+e 
k' 


P 2 X Si{x,k') + x— {x Si{x,k')^ 


= hx 


-k' X 


k' 


d / x^ \ 
^ da; V (1 — a;)^'+^ / 


n(k,£) = 


(1 + /3i + /32)(27ri 

1 


^ j il;{z,y)z^-^/-^dzdy 


(l_^)fc'+2 


(1 — x)^''^'^ dx\ {1 — x) 
[/?2 + x{k' + 1 — /32)] . 


(E3) 


(1 + Pi +/32)(27rf)2 
1 


1 


^k-i3i-iye-/32-i 


- ^dzdy Replacing x with .. , a and k' by fc —di and insert- 

1 — qiz ^ — q2y ^ t- b i+g i+fe 

k-p^ 1 -P 2 result into (Ell, we get 

^ ^ _ In _ ^—dzdv 

(1+/3i+^2)(27ri)2 J J zy-yqi-zq2 1 r P 2 (h R \ i q W 

" TWft - *'l 

(l + /3i+/32)(27ri)2 J J y-q2 (1 + + 3 \ l 3 \\ 

! MI \ ^ + [A + TTft+A- A + > - A>1 

- I dzdv ti I a I a ^/c-/3l+2 Q n„ 1 

(E4) 


1 


(1 +/3i +/32)(27ri) J y - q 2 \y - q 2 J = (1 +/3i +/32)'|^ [ ^(^ + 2) 


k—^1 

9i 


(1+/3i+/32)(27ri) / {y - q 2 Y 


y 


k—Pi+t—f. 


(l + /3^)fc-/3i+2 + 


-dzdy 


Plugging this into (El), we get 


k—j3\ 

dl 


1 {k — Pi+£—P2)\ 1 -P 2 

'^2 


{I + P1 + P2) {k-Pi)\ {£-P2)\ 


(1 + Pi + P 2 ) V k — Pi 


(D5) 


g(k) = + 

^ ^ 1 + Pi 


(E5) 
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